Effective field theory analysis of 3D random field Ising model on isometric lattices 
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Ising model with quenched random magnetic fields is examined for single Gaussian, bimodal and 
double Gaussian random field distributions by introducing an effective field approximation that takes 
into account the correlations between different spins that emerge when expanding the identities. 
Random field distribution shape dependencies of the phase diagrams and magnetization curves are 
investigated for simple cubic, body centered and face centered cubic lattices. The conditions for the 
occurrence of reentrant behavior and tricritical points on the system are also discussed in detail. 
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I. INTRODUCTION 



i Ising model [H, (H which was originally introduced as a model describing the phase transition properties of ferromag- 
' netic materials has been widely examined in statistical mechanics and condensed matter physics. In the course of time, 
basic conclusions of this simple model has been improved by introducing new concepts such as disorder effects on the 
' critical behavior of the systems in question. Ising model in a quenched random field (RFIM) which has been studied 
, over three decades is an example of this situation. The model which is actually based on the local fields acting on the 
. . ■ lattice sites which are taken to be random according to a given probability distribution was introduced for the first 
^ ' time by Larkin 3 for superconductors and later generalized by Imry and Ma [4] . Lower critical dimension dc of RFIM 
has been remained as an unsolved mystery for many years and various theoretical methods have been introduced. 
^ I For example, domain wall argument of Imry and Ma [4| suggests that a transition should exist in three and higher 
■ dimensions for finite temperature and randomness which means that dc = 2 5-8]. On the contrary, dimensional 
reduction arguments Q conclude that the system should not have a phase transition at finite temperature in three 
dimensions or fewer, so dc = 3 [lol-[l3l|. On the other hand, Frontera and Vives [31 showed that a two dimensional 
ferromagnetic RFIM with Gaussian random field distribution exhibits order at zero temperature. 

At the same time, a great many of theoretical and experimental works have paid attention on the RFIM and 
quite noteworth y re sults have been obtained. For instance, it has been shown that diluted antiferromagnets such as 
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Fea;Zni_a;F2[ll, 111], Rb2COa;Mgj_^F4[17|, ll8[ and Co2,Zni_2:F2[18| in a uniform magnetic field just correspond to a 
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ferromagnet in a random uniaxial magnetic field [19l . |20| . Following studies have been devoted to investigate the phase 
diagrams of these systems in depth and in the mean field level it was found that different random field distributions lead 
to different phase diagrams for infinite dimensional models. For example, using a Gaussian probability distribution 
Schneider and Pytte 21] have shown that phase diagrams of the model exhibit only second order phase transition 
properties. On the other hand, Aharony j2^] and Matthis [23i] have introduced bimodal and trimodal distributions, 
respectively and they have reported the observation of tricritical behavior. In a recent series of papers, phase transition 
properties of infinite dimensional RFIM with symmetric double [23] and triple |25| Gaussian random fields have also 
been studied by means of a replica method and a rich variety of phase diagrams have been presented. The situation 
has also been handled on three dimensional lattices with nearest neighbor interactions by a variety of theoretical works 
such as effective field theory (EFT) [H^ll, Monte Carlo (MC) simulations [Ml-lli], pair approximation (PA) [HI and 
series expansion (SE) method [11]. By using EFT, Borges and Silva [1^ studied the system for square (g = 4) and 
simple cubic {q = 6) lattices and they observed a tricritical point only for (g > 6). Similarly, Sarmento and Kaneyoshi 
[2^ investigated the phase diagrams of RFIM by means of EFT with correlations for a bimodal field distribution and 
they concluded that reentrant behavior of second order is possible for the system with (g > 6). Recently, Fytas et al. 
[33I I applied Monte Carlo simulations on a simple cubic lattice. They found that the transition remains continuous 
for a bimodal field distribution, while Hadjiagapiou |42i] observed reentrant behavior and confirmed the existence of 
tricritical point for an asymmetric bimodal probability distribution within the mean-field approximation based on 
Landau expansion. 

Conventional EFT approximations include spin-spin correlations resulting from the usage of the Van der Waer- 
den identities and provide results that are superior to those obtained within the traditional MET. However, these 
conventional EFT approximations are not sufficient enough to improve the results, due to the usage of a decoupling 
approximation (DA) that neglects the correlations between different spins that emerge when expanding the identities. 
Therefore, taking these correlations into consideration will improve the results of conventional EFT approximations. 
In order to overcome this point, recently we proposed an approximation that takes into account the correlations be- 
tween different spins in the cluster of considered lattice |37h41| ] . Namely, an advantage of the approximation method 
proposed by these studies is that no decoupling procedure is used for the higher order correlation functions. On the 
other hand, as far as wc know, EFT studies in the literature dealing with RFIM are based only on discrete probability 
distributions (bimodal or trimodal). Hence, in this work we intent to study the phase diagrams of the RFIM with 
single Gaussian, bimodal and double Gaussian random field distributions on isometric lattices. 

Organization of the paper is as follows: In section|lT]we briefly present the formulations. The results and discussions 
are presented in section [TTTl and finally section ITVl contains our conclusions. 



In this section, we give the formulation of present study for simple cubic lattice with q — 6. A detailed explanation 
of the method for body centered and face centered cubic lattices can be found in Appendix. As a simple cubic lattice, 
we consider a three-dimensional lattice which has N identical spins arranged. We define a cluster on the lattice which 
consists of a central spin labeled Sq, and q perimeter spins being the nearest neighbors of the central spin. The cluster 
consists of (q + 1) spins being independent from the spin operator S. The nearest neighbor spins are in an effective 
field produced by the outer spins, which can be determined by the condition that the thermal average of the central 
spin is equal to that of its nearest neighbor spins. The Hamiltonian describing our model is 



where the first term is a summation over the nearest neighbor spins with 5f — ±1 and the second term represents 
the Zeeman interactions on the lattice. Random magnetic fields are distributed according to a given probability 
distribution function. Present study deals with three kinds of field distribution. Namely, a normal distribution which 
is defined as 



II. FORMULATION 




(1) 



<i,j> 




(2) 



with a width a and zero mean, a bimodal discrete distribution 



(3) 
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where half of the lattice sites subject to a magnetic field ho and the remaining lattice sites have a field —Hq, and a 
double peaked Gaussian distribution 



1 



27rCT" 



1/2 



exp 



ih^-ho)' 

2a2 



- exp 



(4) 



In a double peaked distribution defined in equation random fields ±/io are distributed with equal probability and 
the form of the distribution depends on the Hq and a parameters, where a is the width of the distribution. 

According to Callen identity [52] for the spin-1/2 Ising ferromagnetic system, the thermal average of the spin 
variables at the site i is given by 



{/ijtanh 



(5) 



where j expresses the nearest neighbor sites of the central spin and {/;} can be any function of the Ising variables as 
long as it is not a function of the site. From equation ([5]) with fi — 1, the thermal and random-configurational averages 
of a central spin can be represented on a simple cubic lattice by introducing the differential operator technique [ssl Im ] 



/q=6 



mo = {{S^))r 



l[ [cosh(JV) + 51 sinh(JV)] F{x)\,^o, 



(6) 



where V is a differential operator, q is the coordination number of the lattice and the inner (...) and the outer 
brackets represent the thermal and configurational averages, respectively. The function F{x) in equation ^ is defined 

by 



dhiP{hi) tanh[/3(a; + hi)], 



(7) 



and it has been calculated by numerical integration and by using the distribution functions defined in equations ([2]), 
and Q. By expanding the right hand side of equation ^ we get the longitudinal spin correlation as 



{{So))r = fco + 6fci((5i)), + 15fc2((5i52)). + 20fc3((5l5253))r + lhk^{{SiS2S:,S^)) , 
+&k^{{SiS2S^SAS<s))r + fc6(('5'lS'25354S'5S'6))r- 



(8) 



The coefficients in equation ([SJ are defined as follows 

cosh^( JV) sinh( JV)F(x)U=o 
cosh''(JV)sinh^(JV)F(x) 



fco 
fci 
k2 



x=0, 



ke 



cosh^(JV) sinh"'(JV)F(a;)U=o, 
cosh2(JV) sinh''(JV)F(a;)U=o, 
cosh(JV) sinh^( JV)i^(x)U=o, 
sinhS(JV)F(a;)U=o- 

Next, the average value of the perimeter spin in the system can be written as follows and it is found as 

mi = {{Si))r = ((cosh( JV) + 5*0 sinh( JV)}}^F(x + 7), 

= fli + a2((5o))r- 

For the sake of simplicity, the superscript z is omitted from the left and right hand sides of equations ([8]) and (flO 
The coefficients in equation ([TU]) are defined as 



(9) 



(10) 



fli = cosh(JV)F(a; + 7)|:^:=o, 
02 = sinh(JV)F(a; + 7)|x=o- 



(11) 



In equation (|lll) . 7 = (? — 1)^ is the effective field produced by the {q — 1) spins outside of the system and A is 
an unknown parameter to be determined self-consistently. Equations ([5]) and (jlO[) are the fundamental correlation 
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functions of the system. When the right-hand side of equation ([6]) is expanded, the multispin correlation functions 
appear. The simplest approximation, and one of the most frequently adopted is to decouple these correlations 
according to 

{{S.S,...S^))^^{{S,))^{{S,))^...m),, (12) 

for i ^ j ... I |55| . The main difference of the method used in this study from the other approximations in the 
literature emerges in comparison with any decoupling approximation (DA) when expanding the right-hand side of 
equation (jG]). In other words, one advantage of the approximation method used in this study is that such a kind of 
decoupling procedure is not used for the higher order correlation functions. For spin-1/2 Ising system in a random 
field, taking equations ([8]) and ([TO)) as a basis we derive a set of linear equations of the spin correlation functions in 
the system. At this point, we assume that (z) the correlations depend only on the distance between the spins and 
(ii) the average values of a central spin and its nearest-neighbor spin (it is labeled as the perimeter spin) are equal to 
each other with the fact that, in the matrix representations of spin operator S, the spin-1/2 system has the property 
(S)'^ = 1 . Thus, the number of linear equations obtained for a simple cubic lattice (q = 6) reduces to twelve and the 
complete set is as follows 



{{Sl))r 


= ai+ a2{{So))r, 




= ai{{Si))r + a2{{SoSi))r, 




= ai{{SiS2))r + a2{{SoSiS2))r, 




= ai{{SiS2S3))r + a2{{SoSiS2S3))r, 


{{SiS2S3SiS5)), 


— ai{{SiS2S3S4))r + a2(('S'oS'lS'2<S'3S'4))r, 


( (<5'i 5*2 5*3 5*4 5*5 S'e) ) , 


= ai{{SlS2S3S4S5))r + a2{{SoSiS2S3SiS5))r, 


{{So))r 


= ko+ 6h{{Si))r + 15k2{{SiS2))r + 20A:3 ( (5i ^s^s)) r + l5h({SiS2S3Si))r 




+ 6k5{{SiS2S3S4:S5))r + ke{{SiS2S3SiS5S6))r, 




- 6fcl + (fco + 15fc2)((5l)),. + 20k3{{SlS2))r + 15ki{{SlS2S3))r 




+6k5{{SiS2S3S4))r + ke{{SiS2S3S4,S5))r, 




= (6fci + 20k3){{Si))r + (fco + 15fc2 + 15fc4)((S'l52))r 




+ 6fc5((S'iS'2S'3))r + k6{{SiS2S3S4))r, 




- (6fci + 20fc3 + 6fc5)((5i52))r + (fco + 15fc2 + 15fc4 + h){{S,S2S3))r, 




= (6fci + 20fc3 + 6k5){{SiS2S3))r + (fc-0 + 15fc2 + 15fc4 + h){{SiS2S3Si))r, 


{{SoSiS2S3S4S^))^ 


= (6fci + 20fc3 + 6k5){{SiS2S3S4))r + (fco + 15fc2 + 15fc4 + ke){{SiS2S3S4S5))r 



If equation (|13p is written in the form of a 12 x 12 matrix and solved in terms of the variables Xi[{i — 
1,2, ...,12){e.g.,xi = {{Si))r,X2 = {{SiS2))r,X3 — {{SiS2S3))r, ■■■)] of the linear equations, aU of the spin corre- 
lation functions can be easily determined as functions of the temperature and Hamiltonian/random field parameters. 
Since the thermal and configurational average of the central spin is equal to that of its nearest-neighbor spins within 
the present method, the unknown parameter A can be numerically determined by the relation 

{{So))r = {{Si))r or X7^Xi. (14) 

By solving equation (fT4|) numerically at a given fixed set of Hamiltonian/random field parameters we obtain the 
parameter A. Then we use the numerical values of A to obtain the spin correlation functions which can be found 
from equation (|13|) . Note that A = is always the root of equation (IT4|) corresponding to the disordered state of 
the system. The nonzero root of A in equation ()14|) corresponds to the long range ordered state of the system. 
Once the spin correlation functions have been evaluated then we can give the numerical results for the thermal and 
magnetic properties of the system. Since the effective field 7 is very small in the vicinity of fcsTc/ J, we can obtain the 
critical temperature for the fixed set of Hamiltonian/random field parameters by solving equation (jl4[) in the limit of 
7—^0 then we can construct the whole phase diagrams of the system. Depending on the Hamiltonian/random field 
parameters, there may be two solutions (i.e. two critical temperature values satisfy the equation (|14p ') corresponding 
to the first/second and second order phase transition points, respectively. We determine the type of the transition by 
looking at the temperature dependence of magnetization for selected values of system parameters. 
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III. RESULTS AND DISCUSSION 

In this section, we discuss how the type of random field distribution effects the phase diagrams of the system. 
Also, in order to clarify the type of the transitions in the system, we give the temperature dependence of the order 
parameter. 

A. Phase diagrams of single Gaussian distribution 

The form of single Gaussian distribution which is defined in equation ^ is governed by only one parameter a 
which is the width of the distribution. In Fig. 1, we show the phase diagram of the system for simple cubic, body 
centered cubic and face centered cubic lattices in {kBTc/J — cr) plane. We can clearly see that as a increases then 
the width of the distribution function gets wider and randomness effect of magnetic field distribution on the system 
becomes significantly important. Therefore, increasing a value causes a decline in the critical temperature ksTc/ J 
of the system. We note that the critical temperature of the system reaches to zero at a = 3.8501, 5.450 and 8.601 
for simple cubic {q = 6), body centered (q = 8) and face centered cubic {q = 12) lattices, respectively. Besides, we 
have not observed any reentrant /tricritical behavior for single Gaussian distribution, or in other words, the system 
undergoes only a second order phase transition. ksTc/J value in the absence of any randomness i.e. when cr = is 
obtained as ksTc/J = 4.5274, 6.5157 and 10.4986 for q = 6,8,12, respectively. These values can be compared with 
the other works in the literature. In Table 1, we see that present work improves the results of the other works. Hence, 
it is obvious that the present method is superior to the other EFT methods in the literature. The reason is due to 
the fact that, in contrast to the previously published works mentioned above, there is no uncontrolled decoupling 
procedure used for the higher order correlation functions within the present approximation. 



TABLE I: Critical temperature ksTc/ J at ho/ J — and cr = obtained by several methods and the present work 
for (? = 6, 8, 12. 



Lattice EBPAf43] 


CEFTf44] 


PA [45] 


EFTf46] 


BAf47] 


EFRGf48] MFRG[49] MC[50] 


SEf51] Present Work 


SC 4.8108 


4.9326 


4.9328 


5.0732 


4.6097 


4.85 4.93 4.51 


4.5103 4.5274 


BCC 


6.9521 








6.88 6.95 6.36 


6.3508 6.5157 


FCC 


10.9696 










9.7944 10.4986 



Magnetization surfaces and their projections on (fc^T/J — cr) plane corresponding to the phase diagrams shown 
in Fig. 1 are depicted in Fig. 2 with q = 6,8 and 12. We see that as the temperature increases starting from zero, 
magnetization of the system decreases continuously and it falls rapidly to zero at the critical temperature for all a 
values. Moreover, critical temperature of the system decreases and saturation value of magnetization curves remains 
constant for a while then reduces as a value increases. This is a reasonable result since, if a value increases, then the 
randomness effects increasingly play an important role on the system and random fields have a tendency to destruct 
the long range ferromagnetic order on the system, and hence magnetization weakens. These observations are common 
properties of both three lattices. 

B. Phase diagrams of bimodal distribution 

Next, in order to investigate the effect of the bimodal random fields defined in equation ([3]) on the phase diagrams 
of the system, we show the phase diagrams in (fc^Tc/J — ho/J) plane and corresponding magnetization profiles with 
coordination numbers q = 6,8, 12 in Figs. 3 and 4. In these figures solid and dashed lines correspond to second 
and first order phase transition points, respectively and hollow circles in Fig. 3 denote tricritical points. As seen in 
Fig. 3, increasing values of ho/J causes the critical temperature to decrease for a while, then reentrant behavior of 
first order occurs at a specific range of ho/J. According to our calculations, reentrant phenomena and the first order 
phase transitions can be observed within the range of 2.0 <ho/J< 3.0 for q = 6, 3.565 < ho/J< 3.95 for q = 8 
and 4.622 < ho/J < 5.81 for q = 12. If ho/J value is greater than the upper limits of these field ranges, the system 
exhibits no phase transition. The tricritical temperatures ksTt/J which are shown as hollow circles in Fig. 2 are 
found as fc^Tt/ J = 1.5687, 2.4751 and 4.3769 for g = 6, 8, 12, respectively. 

In Fig. 4, we show two typical magnetization profiles of the system. Namely, the system always undergoes a second 
order phase transition for ho/J = 1.0. On the other hand, two successive phase transitions (i.e. a first order transition 



FIG. 1: Phase diagrams of simple cubic (g — 6), body centered cubic (q = 8) and face centered cubic {q = 12) lattices in 
{kBTc/J — (j) plane corresponding to single Gaussian distribution of random fields. The numbers on each curve denote the 
coordination numbers. 




is followed by a second order phase transition) occur at the values of ho/ J — 2.5,3.8 and 5.5 for q — 6,8 and 12, 
respectively, which puts forward the existence of a first order reentrant phenomena on the system. We observe that 
the increasing hg/ J values do not effect the saturation values of magnetization curves. 
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FIG. 3: Phase diagrams of the system (g = 6,8,12) in {ksTcj J — ho /J) plane, corresponding to bimodal random field 
distribution. Solid and dashed lines correspond to second and first order phase transitions, respectively. Hollow circles denote 
the tricritical points. 



(a) (b) (c) 



E 




FIG. 4: Temperature dependence of magnetization corresponding to Fig. 3 with a = Q and ho /J = 1.0,2.5 for simple cubic 
lattice (left panel), ho/J = 1.0,3.8 for body centered cubic lattice (middle panel) and ho/J — 1.0,5.5 for face centered cubic 
lattice (right panel). Solid and dashed lines correspond to second and first order phase transitions, respectively. 

C. Phase diagrams of double Gaussian distribution 

Double Gaussian distribution in equation ^ with nearest neighbor interactions have not yet been examined in 
the hterature. Therefore, it would be interesting to investigate the phase diagrams of the system with random fields 
corresponding to equation (U)). Now the shape of the random fields is governed by two parameters h^j/J and u. 
As shown in preceding results, increasing a value tends to reduce the saturation value of the order parameter and 
destructs the second order phase transitions by decreasing the critical temperature of the system without exposing 
any reentrant phenomena for h^/J = 0. Besides, as ha/ J value increases then the second order phase transition 
temperature decreases again and the system may exhibit a reentrant behavior for u = 0, while saturation value of 
magnetization curves remains unchanged. Hence, the presence of both h^/J and a on the system should produce 
a cooperative effect on the phase diagrams of the system. Fig. 5 shows the phase diagrams of the system in 
{kBTc/J —ho/J) and {ksTc/ J — a) planes for g = 6, 8 and 12. As we can see on the left panels in Fig. 5, the system 
exhibits tricritical points and reentrant phenomena for narrow width of the random field distribution and as the width 
a of the distribution gets wider then the reentrant phenomena and tricritical behavior disappear. In other words, both 
the reentrant behavior and tricritical points disappear as a parameter becomes significantly dominant on the system. 
Our results predict that tricritical points depress to zero at ct = 1.421, 2.238 and 3.985 for g = 6, 8 and 12, respectively. 
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FIG. 5: Phase diagrams of the system for a double Gaussian random field distribution with q = 6,8, 12 in {kBTc/J — ho/ J) 
and (knTc/ J — a) planes. 



For distribution widths greater than these values, all transitions are of second order and as we increase cr value further, 



9 




k T/J k T/J k T/J k T/J 

(a) (b) 



FIG. 6: Magnetization curves for simple cubic lattice corresponding to the phase diagrams in Fig. 5 for a double Gaussian 
distribution with some selected values of ho/ J and a. 

then ferromagnetic region gets narrower. Similarly, on the right panels in Fig. 5, we investigate the phase diagrams 
of the system in (fc^Tc/J — a) plane with selected values of ho/ J. We note that for the values oi ho/ J < 2.0 {q — 6), 
ho/J< 3.565 {q ~ 8) and ho/ J < 4.622 {q = 12), system always undergoes a second order phase transition between 
paramagnetic and ferromagnetic phases at a critical temperature which decreases with increasing values of ho/ J as 
in Fig. 1, where ho/ J = 0. Moreover, for the values of ho/ J greater than these threshold values the system exhibits 
a reentrant behavior of first order and the transition lines exhibit a bulge which gets smaller with increasing values 
of ho/J which again means that ferromagnetic phase region gets narrower. Besides, for ho/ J > 2.9952 (g = 6), 
ho/J > 3.9441 {q — 8) and ho/J > 5.8085 {q — 12) tricritical points appear on the system. In Fig. 6, we show 
magnetization curves corresponding to the phase diagrams shown in Fig. 5 for simple cubic lattice. Fig. 6a shows 
the temperature dependence of magnetization curves for q — 6 with ho/J — 0.5 (left panel) and ho/J = 2.5 (right 
panel) with selected values of a. As we can see in Fig. 6a, as a increases then the critical temperature of the system 
decreases and first order phase transitions disappear (see the right panel in Fig. 6a). Moreover, rate of decrease of the 
saturation value of magnetization curves increases as ho/J increases. On the other hand, in Fig. 6b, magnetization 
versus temperature curves have been plotted with cr = 0.5 (left panel) and cr = 2.5 (right panel) with some selected 
values of ho/J. In this figure, it is clear that saturation values of magnetization curves remain unchanged for a = 0.5 
and tend to decrease rapidly to zero with increasing values of ho/J when a — 2.5. In addition, as ho/J increases 
when the value of a is fixed then the critical temperature of the system decreases and ferromagnetic phase region of 
the system gets narrower. These observations show that there is a collective effect originating from the presence of 
both ho/J and a parameters on the phase diagrams and magnetization curves of the system. 

Finally, Fig. 7 represents the variation of the tricritical point {ksTt/ J, ht/ J) with a for q = 6,8 and 12. As seen 
from Fig. 7, ksTt/ J value decreases monotonically as a increases and reaches to zero at a critical value at- According 
to our calculations, the critical distribution width at value can be given as at = 1.421,2.238,3.985 for q — 6,8 and 
12, respectively. This implies that at value depends on the coordination number of the lattice. Besides, ht/ J curves 
exhibit relatively small increment with increasing a value. 

IV. CONCLUSIONS 

In this work, we have studied the phase diagrams of a spin- 1/2 Ising model in a random magnetic field on simple 
cubic, body centered cubic and face centered cubic lattices. RFIM with nearest neighbor interactions with double 
Gaussian distribution as well as its single Gaussian counterpart has not been examined within the EFT approximation 
before. Hence, we have introduced an effective field approximation that takes into account the correlations between 
different spins in the cluster of considered lattice and examined the phase diagrams as well as magnetization curves 
of the system for different types of random field distributions, namely single Gaussian, bimodal and double Gaussian 
distributions. For a single Gaussian distribution, we have found that the system always undergoes a second order phase 
transition between paramagnetic and ferromagnetic phases. In the absence of any randomness (i.e. ho/ J — 0,a — 0) 
critical temperature values corresponding to the coordination numbers q = 6,8 and 12 are found to be in an excellent 
agreement with those of the other theoretical methods such as MC and SE. For bimodal and double Gaussian 
distributions, we have given the proper phase diagrams, especially the first-order transition lines that include reentrant 
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phase transition regions. Our numerical analysis clearly indicate that such field distributions lead to a tricritical 
behavior. Moreover, we have discussed a collective effect which arise from the presence of both ho/ J and a parameters 
together and we have observed that saturation values of magnetization curves are strongly related to this collective 
effect. 

As a result, we can conclude that all of the points mentioned above show that our method is superior to the 
conventional EFT methods based on decoupling approximation and therefore we think that the calculation results 
are more accurate. Therefore, we hope that the results obtained in this work may be beneficial from both theoretical 
and experimental points of view. 
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Appendix 

The number of distinct correlations for the spin- 1/2 system with q nearest neighbors is 2g — 1. With the mag- 
netization of the central spin {{So))r and its nearest neighbor spin {{Si))r, there has to be 2g -I- 1 unknowns which 
forms a system of linear equations, q of the correlations which include central spin are derived from the central spin 
magnetization {{So))r, and rest of the q — 1 correlations which include only perimeter spin are derived from perimeter 
spin magnetization expression {{Si))r- Let us label these correlations with Xi, i = 1, 2, . . . , 2g -I- 1, such that first q of 
them include only perimeter spins and last q + 1 of them include central spin. We can represent all of the spin corre- 
lations and centeral/perimeter spin magnetization with Xi as xi = {{Si))r, X2 = {{SiS2))r,- ■ -Xq = ((5*1 5*2 . . . Sq))r, 
Xq+i = {{So))r, Xq+2 = {{SoSi))r, ■ ■ ■ and X2q+i — {{SqSi . . . Sq))r- The Complete set of correlation functions for 
q — 8, 12 can be obtained as the same procedure given between equations (p|)- (IT51) in section HIl 

The complete set of correlation functions for body centered cubic lattice {q — 8): 

xi = ai + a2Xg, 

X2 = aixi + 02X10, 

X3 = aiX2 + 02X11, 

X4 = 01X3 -I- 02X12, 
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X5 = aiXi + a2Xi3, 

xq = aiX5 + a2Xi4, 

X7 = aixe + a2Xi5, 

xs = aiXr + a2Xie, 

xg = ko + SkiXi + 28^2X2 + 56^30:3 + 70/;;4a;4 + bGk^x^ + 28kQXQ + SkjXr + kgxs, 

xio = 8fci + (28fc2 + kf))xi + 56fc3a;2 + 70fc4.X3 + 56^:5X4 + 28^6X5 + SkjXa + k^XT, 

Xii = (56^3 + 8fci)xi + (70fc4 + 28fc2 + A;o)x2 + 56fc5X3 + 28/05X4 + 8/07X5 + fcgXe, 

xi2 = (56^3 + 56fc5 + 8fci)x2 + (70fc4 + 28^2 + A;o + 28fc6)x3 + 8^7x4 + ^8X5, 

xi3 = (56fc3 + 56fc5 + 8fc7 + 8A;i)x3 + (70fc4 + /;;8 + 28A;2 + fco + 28fc6)a;4, 

xi4 = (56^3 + 56^5 + 8^7 + 8A;i)x4 + (70fc4 + A;8 + 28A;2 + fco + 28fc6)a;5, 

xi5 = (56fc3 + 56fc5 + 8^7 + 8A;i)x5 + (70^4 + fcs + 28fc2 + fco + 28A:6)a;6, 

X16 = (56^3 + 56fc5 + 8fc7 + 8A;i)x6 + (70fc4 + fcs + 28fc2 + ko - 28ke)x7, 

xi7 = (56^3 + 56^5 + 8^7 + 8A;i)x7 + (70fc4 + A:8 + 28/02 + fco + 28/c6)x8. 



Similarly, correlation functions for face centered cubic lattice {q = 12): 

xi = ai+ 02X13, 

X2 = aiXi + a2Xi4, 

X3 = aiX2 + 02X15, 

X4 = 01X3 + 02X16, 

X5 = aiX4 + 02X17, 

X6 = 01X5 + 02X18, 

X7 = OiXe + 02X19, 

X8 = 01X7 + 02X20, 

Xg = aiX8 + 02X21, 
Xio = OiXg + O2X22, 
Xii = OiXio + 02X23, 
X12 = OiXii + 02X24, 

xi3 = /oo + 12/eiXi + 66/02X2 + 220/03X3 + 495/04X4 + 792/05X5 + 924/06X6 

+792/07X7 + 495/08X8 + 22O/09X9 + 66/oioXio + 12fciiXii + /012X12, 
xi4 = 12A;i + (fco + 66fc2)xi +220^3X2 +495/04X3 + 792A;5X4 + 924/062:5, 

+792/07X6 + 495/08X7 + 22O/09X8 + 66/oioa;9 + 12/oiiXio + ki2Xu, 
xi5 = (12fci + 220/O3)xi + (fco + 66/02 + 495fc4)x2 + 792/05X3 + 924^6X4, 

+792/07X5 + 495/08X6 + 22O/09X7 + 66/010X8 + I2/011X9 + /012X10, 
X16 = (12fci+ 220/03 + 792/o5)x2 + (/oo + 66/02 + 495A;4 + 924/o6)x3, 

+792A;7X4 + 495/o8X5 + 220^9X5 + GGkioXr + 12fciiX8 + A;i2a;9, 
xi7 = (12fci + 220^3 + 792/05 + 792fc7)x3 + {ko + 66^2 + 495^4 + 495/o8 + 924A;6)x4, 

+22O/09X5 + 66/010X6 + I2/011X7 + /012X8, 
X18 = (12fci+ 220/03 + 220/09 + 792/05 + 792/o7)x4, 

+(/oo + 66/02 + 495/04 + 495/08 + 66fcio + 924fc6)x5 + 12/oiiX6 + /oi2a;7, 
xi9 = (12fci+ 220/03 + 220/09 + 792^5 + 12fcii + 792fc7)x5, 

+(/oo + 66/02 + A;i2 + 495/04 + 495/08 + 66A;io + 924fc6)x6, 
X20 = (12fci+ 220/03 + 220/09 + 792/05 + 12fcii + 792/o7)x6, 

+(/oo + 66/02 + /012 + 495/04 + 495/08 + 66/010 + 924/o6)x7, 
X21 = (12fci + 220/03 + 220/09 + 792^5 + 12fcii + 792fc7)x7, 
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+ (fco + 66fc2 + fci2 + 495^4 + 495fc8 + 66fcio + 924fc6)x8, 
X22 = (12fci+220fc3 + 220fc9 + 792^-5 + 12A;ii+ 792fc7)a;8, 

+(A;o + 66A;2 + ki2 + 495A;4 + 495fc8 + 66fcio + 924fc6)x9, 
X23 = (12fci+220A;3 + 220A;9 + 792^5 + 12fcii + 792^7)0:9, 

+{ko + 66k2 + ki2 + 495A;4 + 495A;8 + 66fcio + 924fc6)a;io, 
a;24 = (12fci +220^3 + 220^9 + 792^5 + 12fcii + 792fc7)a;io, 

+(fco + 66fc2 + ki2 + 495^4 + 495fc8 + 66fcio + 924fc6)a;ii, 
X25 = (12fci+ 220/03 + 220/09 + 792^5 + 12fcii + 792fc7)a;ii, 

+{ko + 66A;2 + ki2 + 495/04 + 495^8 + 66fcio + 924fc6)a;i2. 

where the coefficients are given by 

ao = cosh.{JV)F{x + 7)|x=o, cti = smk{JV)F{x + 7)|a;=o, 

with 7 = (g - 1)^, 

ki = cosh«-'(JV) sinh'(JV)F(a;)|x=o, i = 0,...,q (15) 

where q is the coordination number. Since the first 2q equations are independent of the correlation labeled X2q+i it 
is not necessary to include this correlation function in the set of linear equations. Therefore it is adequate to take 2q 
hnear equations for calculations. 
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